library(raster)
library(ncdf)
library(gtools)



## specify outpath for storage of renamed files
dirpath <- "E:/TRMMDaten/incoming/TRMM04bis06"
outpath <- "E:/TRMMDaten/incoming/TRMM04bis06Out"

station <- read.csv("E:/Messstationsdaten/add_3hour/april/ki_0000cof3_000pu1_201104010000_201104302355_eat_ca01_nai05_0010.dat_suh03")

dirpath_stat <- "E:/Messstationsdaten/add_3hour/cof3"
  
## list all files in dirpath with that have '.nc' in their name
fnames <- list.files(dirpath, recursive = T, pattern = ".nc")

## extract long and short file names
##nchar zaehlt die Anzahl der Elemente des Vektors
fnames_long <- subset(fnames, nchar(fnames) == 20)
fnames_short <- subset(fnames, nchar(fnames) == 19)

## extract hour from short names and create long hour specifications
##substr gibt einen Teil des Vektors zur??ck
##paste verbindet Vektoren nachdem sie in characters umgewandelt worden sind
short_hour <- substr(fnames_short, 13, 13)
long_hour <- paste("0", short_hour, sep = "")
short_hour <- paste(".", short_hour, ".", sep = "") ##wofuer?
long_hour <- paste(".", long_hour, ".", sep = "")

## replace hour bit in short file names with long hour specification
##gsub ersetzt ein Textmuster durch eine Zeichenfolge im GESAMTEN Text
fnames_short <- lapply(1:length(fnames_short), function(i) {
  gsub(short_hour[i], long_hour[i], fnames_short[i], fixed = T)
}
                            )
## create new fnames_new list (with long hour specification)
fnames_new <- c(fnames_long, fnames_short)
fnames_new <- mixedsort(fnames_new) ## sort into proper temporal order
##mixedsort sortiert Strings in die richtige Reihenfolge

## now rename files on disk so they comply with new (longer) file names
file.copy(from = paste(dirpath, fnames, sep = "/"), 
          to = paste(outpath, fnames_new, sep = "/"))

## list all station files
fnames_stat <- list.files(dirpath_stat, recursive = T, pattern = "ki")

## read in all files in 'names_list' - returns list with data frames
data_list <- lapply(1:length(fnames_new), function(i) {
  open.ncdf(paste(outpath, fnames_new[i], sep = "/"), write = F, readumlim = T)
  }
                    )
data_list_stat <- lapply(1:length(fnames_stat), function(i) {
  read.table(paste(dirpath, fnames[i], sep = "/"), header = F, sep = ",", fill = T, skip = 1)
}
                         )


##get variable Precipitation and transpose columns and rows
Prec <- lapply(1:length(data_list), function(i) {
  t(get.var.ncdf(data_list[[i]], varid = "precipitation"))
  }
               )

##create Raster

PrecRaster <- lapply(1:length(Prec), function(i) {
  raster(Prec[[i]], xmn = -180, xmx = 180, ymn= -50, ymx = 50)
})


##extract values of raster list (no raster any more)
PrecExtract <- lapply(1:length(PrecRaster), function(i) {
    extract(PrecRaster[[i]], extent(37, 37.25, -3.25, -3))
}
              )

##to get the correct amount of rainfall multiply the values by 3 
##(TRMM Data in mm h^-1)
PrecExtractProd <- lapply(1:length(PrecExtract), function(i) {
  (PrecExtract[[i]] * 3)
})

##plot 
xyplot(unlist(PrecExtractProd) ~ unlist(data_list_stat$precip_pu1), type="p", 
       panel = function(x, y, ...){
  panel.xyplot(x, y, ...)
  panel.lmlineq(x, y, rot=T, at=0.7, r.squared = T)
})
